package scu.maqiang.cpte;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.fes.FES1L2CTE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

public class CoupledPiezoThermalElasticity {

    public static void main(String[] args) {
        double kxx = 1.0, density = 1.0, specificHeat = 1.0;
        double E = 2.0, Nu = 0.3, alpha = 0.05;
        double beta = E / (1 - 2 * Nu)  * alpha;

        double theta = 0.5, T = 1.5e1, dt = T / 1500;
        double omega = 0.25;

        Mesh1L2 mesh2D = new Mesh1L2();
        mesh2D.line(0, 1.0, 100);

        FES1L21 ffs = new FES1L21(mesh2D);

        FES1L2CTE fs = new FES1L2CTE(mesh2D, ffs);

        SRMatrix L = new SRMatrix(fs.GetNdof());
        SRMatrix R = new SRMatrix(fs.GetNdof());
        SRMatrix ElasticMass = new SRMatrix(ffs.GetNdof());
        SRMatrix HTEMatrix = new SRMatrix(ffs.GetNdof(), ffs.GetNdof());

        double[] x = new double[fs.GetNdof()];
        double[] RHS = x.clone();

        double[] temp = new double[ffs.GetNdof()];
        double[][] dispUV = new double[2][ffs.GetNdof()];

        double[] totalDisp = new double[ffs.GetNdof()];
        double[] totalDisp0 = totalDisp.clone();
        double[] vel = totalDisp.clone();
        double[] vel0 = totalDisp.clone();
        double[] acc = totalDisp.clone();
        double[] acc0 = totalDisp.clone();

        double[] ex = temp.clone();
        double[] ey = temp.clone();
        double[] exy = temp.clone();
        double[] ez = temp.clone();
        double[] tx = temp.clone();
        double[] ty = temp.clone();
        double[] txy = temp.clone();

        double MM = E * (1 - Nu) / ((1 + Nu) * (1 - 2 * Nu));
        double TRef = 8.0 * 35.0 / 13.0;
        double AA = MM + TRef * beta * beta / (density * specificHeat);
        System.out.println("A = " + AA);
        System.out.println("M = " + MM);

        ScalarFunc initTemp = (xy, label, param) -> xy[0] > 0.999999? 1.0:0.0;
        double[] Itemp = ffs.valueFromFunc(initTemp, null);
        ffs.assembleMass(new double[]{density}, BVPType.CONSIST_MASS, ElasticMass);
        fs.assembleHeatToElastic(beta, BVPType.COMMON, HTEMatrix);
        HTEMatrix.mul(Itemp, totalDisp);
        NewIterSSolver solver = new NewIterSSolver(ElasticMass);
        solver.PCGSSOR(totalDisp, acc, 1.5, 1);
        Arrays.setAll(acc0, i -> acc[i]);
        double[] coef = {density, specificHeat, kxx, TRef, theta, E, Nu, beta, dt, omega};
        fs.assembleDierctCouplingMatrix_Type2(coef, BVPType.COMMON, BVPType.CONSIST_MASS,
                BVPType.COMMON, BVPType.CONSIST_MASS, BVPType.COMMON, L, R);
        //L.sort();
        //System.out.println(L);

        fs.applyBCHeat_MBN(L, 2);
        fs.applyBCDisp_MBN(L, Direct.All, 1);
        //System.out.println(L);
        solver.setMatrix(L);
        double time = 0.0;
        double[] RHS1 = new double[fs.GetNdof()];

        int lineNum = 100;
        double dis = 1.0 / lineNum;
        double[] lineX = new double[lineNum + 1];
        double[] lineY = lineX.clone();

        Arrays.setAll(lineX, i -> i * dis);
        String dir = "RectangleResultL2\\";
        String XYLineFile = dir + "Line150_2.dat";
        Tecplot.lineXYTime(XYLineFile, lineX, lineY, 0.0);

        String heatFile = dir + "Heat150_2.dat";
        mesh2D.toTecplot(heatFile, 0.0, x);

        String deformFile =  dir + "Deform150_2.dat";
        mesh2D.toTecplotDeform(deformFile, 0.0, dispUV, 0.0);

        String expansionFile = dir + "expansion150_2.dat";
        String temp00File = dir + "temp00_150_2.dat";
        FileWriter fw1 = null;
        FileWriter fw2 = null;
        try {
            fw1 = new FileWriter(expansionFile);
            fw2 = new FileWriter(temp00File);
            fw1.write("0.0\t0.0\n");
            fw2.write("0.0\t0.0\n");
        } catch (IOException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
        while(time < T) {
            time += dt;
            System.out.println("time: " + time);
            fs.applyBCHeat_MBN(RHS, 1.0, 2);
            fs.applyBCDisp_MBN(RHS, Direct.All, 0.0, 1);
            solver.PCGSSSOR(RHS, x, 1.5, 1);
            fs.extractUVH(x, temp, dispUV);
            mesh2D.toTecplot(heatFile, time, temp);
            System.out.println("T min: " + MVO.min(temp)  + "\t max: " + MVO.max(temp));
            System.out.println("U min: " + MVO.min(dispUV[0]) + "\t max: " + MVO.max(dispUV[0]));

            for(int i = 0, n = ffs.GetNdof(); i < n; i++) {
                totalDisp[i] = dispUV[0][i];
            }

            Arrays.setAll(acc, i -> 1/(omega*dt*dt) * (totalDisp[i] - totalDisp0[i]) -1/(omega * dt) * vel0[i] - (0.5 /omega - 1) * acc0[i]);
            Arrays.setAll(vel, i -> vel0[i] + ((1-0.5) * acc0[i] + 0.5 * acc[i]) * dt);

            Arrays.setAll(acc0, i -> 1/(omega*dt) * vel[i] + (0.5/omega -1) * acc[i]);
            ElasticMass.mul(acc0, vel0);
            for(int i = 0, n = ffs.GetNdof(); i < n; i++) {
                RHS1[i] = vel0[i];
            }

            R.mul(x, RHS);
            Arrays.setAll(RHS, i -> RHS[i] + RHS1[i]);
            Arrays.setAll(acc0, i -> acc[i]);
            Arrays.setAll(vel0, i -> vel[i]);
            Arrays.setAll(totalDisp0, i -> totalDisp[i]);

            mesh2D.toTecplot(heatFile, time, temp);
            mesh2D.toTecplot(deformFile, time, dispUV[0]);
            Arrays.setAll(dispUV[0], i -> dispUV[0][i] * (1-Nu) / (alpha * (1 + Nu)));
            Tecplot.lineXYTime(XYLineFile, lineX, dispUV[0], time);
            try {
                fw1.write(time + "\t " + dispUV[0][100] * (1-Nu)/(alpha * (1 + Nu))  + "\n");
                fw2.write(time + "\t " + temp[0] + "\n");
            } catch (IOException e) {
                // TODO Auto-generated catch block
                e.printStackTrace();
            }
        }

        try {
            fw1.flush();
            fw1.close();
            fw2.flush();
            fw2.close();
        } catch (IOException e) {
            e.printStackTrace();
        }
    }
}
